clear all
close all
clc

% Output Selection

report       = 'T4';          % 'T1', 'T3', 'T4' 
markup       = '1.15';        % '1.05', '1.15', '1.25', '1.35'
experiment   = 'sizedep';   % 'efficient', 'uniform', 'sizedep', 'entry'

% Options

set(groot, 'DefaultAxesLineWidth', 1.5);
set(groot, 'DefaultLineLineWidth', 3);
set(groot, 'DefaultAxesTickLabelInterpreter','latex'); 
set(groot, 'DefaultLegendInterpreter','latex');
set(groot, 'DefaultAxesFontSize', 22);

optset('bisect','tol', 1e-18); 

options.Display = 'none';

% Calibrated Parameters

switch markup 

    case '1.05'

        Markuptarget = 1.05;                      % our target for aggregate markup
        
        p.xi         = 20.7021016669651;          % Pareto tail productivity
        p.sigma      = 29.1008275310716;          % demand elasticity: match level of markups 

    case '1.15'

        Markuptarget = 1.15;                      % our target for aggregate markup
        
        p.xi         =  6.8409724072258;          % Pareto tail productivity
        p.sigma      = 10.8610835311753;          % demand elasticity: match level of markups 

    case '1.25'

        Markuptarget = 1.25;                      % our target for aggregate markup
        
        p.xi         = 4.06857848423202;          % Pareto tail productivity
        p.sigma      = 7.2079290139342;           % demand elasticity: match level of markups 

    case '1.35'

        Markuptarget = 1.35;                      % our target for aggregate markup
        
        p.xi         = 2.88859963714985;          % Pareto tail productivity
        p.sigma      = 5.66479780771211;          % demand elasticity: match level of markups 

end

% Assigned Parameters

mshare     = 0.45;                        % share of intermediate inputs in total sales

p.epsi     =  0.162*p.sigma; 
p.phi      =  0.4;                        % doesn't matter since reset below to match material share  

p.beta     = 0.96;                        % period is 1 year
p.varphi   = 0.04;                        % exit rate 
p.delta    = 0.06;                        % depreciation rate

p.r        = 1/p.beta - 1; 
p.nu       = 1;                           % inverse labor supply elasticity
p.alpha    = 1/3;                         % capital share in value added
p.theta    = 0.5;                         % elasticity of substitution between VA and intermediates

p.xie      = 0;                           % entry subsidy
p.xis      = 0;                           % sales subsidy
p.taus     = 0; 

% Initial Aggregates (Normalization)

Y          = 1; 
N          = 1; 
R          = p.r + p.delta; 

% Quadrature
 
nz         = 5000;                         % Gauss-Legendre quadrature for exponential r.v.

zmax       = -1/p.xi*log(1e-22);

[z, w]     = qnwlege(nz, 0, zmax); 
w          = p.xi.*w.*exp(-p.xi.*z);
w          = w./sum(w); 

z          = exp(z); 

p.w        = w; 
p.z        = z; 


%%%%% 1. Solve D and firm's choices for a given N %%%%%% 

%%%% Start solving

D              = fsolve('findequilibrium', 1, options, w, z, p, 'market', 'old');

[~, ~, ~, Omega, ~,~, ~, ~, ~, ~, ~, ~, Mu]  = findequilibrium(D, w, z, p, 'market', 'old');

p.phi          = 1 - mshare*Omega^(1 - p.theta)*Mu;                                       % choose to match materials share in sales = mshare 

[~, q, mu, Omega, Z, W, Lp, K, B, U, Uq, C, Mu]  = findequilibrium(D, w, z, p, 'market', 'old');

p.F            = Y/N/W/(1/p.beta - 1 + p.varphi)*(1 - 1/Mu);     % uses normalization N/Y = 1 to back out F

L              = Lp + p.F*p.varphi*N;

p.psi          = W/(C*L^p.nu); 

P              = (N*w'*(Uq.*q))^(-1);       % Demand index so we can compute prices

% Individual Choices

py             = Uq*P;                      % producer price
y              = q*Y; 
l              = py.*y/Y.*Mu./mu.*Lp; 
b              = py.*y/Y.*Mu./mu.*B; 
k              = py.*y/Y.*Mu./mu.*K; 
mc             = Omega./z; 

Dp       = fsolve('findequilibrium', D, options, w, z, p, 'planner', 'old');

[~, qp, ~, ~, Zp, ~,~, ~, ~, Up]  = findequilibrium(Dp, w, z, p, 'planner', 'old');

% compute value added misallocation

GDPp  = p.phi^(1/(p.theta - 1))*Zp/(1 - (1 - p.phi)*Zp^(p.theta - 1))^(1/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
Yp    = p.phi^(1/(p.theta - 1))*Zp/(1 - (1 - p.phi)*Zp^(p.theta - 1))^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
Bp    = (1 - p.phi)*Zp.^(p.theta - 1)*Yp; 

GDP   = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1)*Mu^(-p.theta))*Z/(1 - (1 - p.phi)*Z^(p.theta - 1)*Mu^(1 - p.theta))^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
GDP1  = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1)              )*Z/(1 - (1 - p.phi)*Z^(p.theta - 1)                 )^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);

moment_data        = [Markuptarget; 0.45; 0.229; 0.573; 0.739; 0.162]; 

moment_model       = zeros(size(moment_data)); 

moment_model(1)    = Mu; 
moment_model(2)    = B/Y;

data               = sortrows([w, py.*y], 2); 

data               = data(data(:, 2) > 1e-16, :); 

data(:,1)          = data(:,1)/sum(data(:,1)); 

cumF               = cumsum(data(:,1)); 

top1               = cumF >= 0.99; 
top5               = cumF >= 0.95; 
top10              = cumF >= 0.90; 

moment_model(3)    = data(top1,  1)'* data(top1,  2)/(data(:,1)'*data(:,2));
moment_model(4)    = data(top5,  1)'* data(top5,  2)/(data(:,1)'*data(:,2));
moment_model(5)    = data(top10, 1)'* data(top10, 2)/(data(:,1)'*data(:,2));

keep               = py.*y > 1e-8; 

bhat               = lscov([ones(size(y(keep))), log(py(keep).*y(keep))], 1./mu(keep) + log(1 - 1./mu(keep)));

moment_model(6)    = bhat(2);

weight             = w.*l/sum(w.*l); 

weights            = zeros(numel(moment_model), 1); 

weights([1; 2; 4]) = 1; 

weights            = weights/sum(weights);

err_mom            = (moment_model - moment_data)./(1 + moment_data);
err_mom            = (weights'*err_mom.^2).^(1/2);

profits            = (P*Uq.*q - Omega*(q./z));

transfer           = (U - Uq.*q)*P;

net_profits        = profits+transfer;

% Output 

switch report

    case 'T1'

        symbol      = categorical({sprintf('%c ', 946);sprintf('%c ', 948);sprintf('%c ', 966);sprintf('%c ', 945);sprintf('%c ', 957);sprintf('%c ', 952)});
        description = categorical({'discount factor';'depreciation rate';'exit rate';'elasticity of value-added to capital';'elasticity of labor supply';'elasticity of substitution between value-added and materials'});
        values      = [p.beta; p.delta; p.varphi; p.alpha; p.nu; p.theta];
        
        fprintf('<strong> Table 1: Parameterization, Benchmark Model </strong> \n'); 

        fprintf('\n');
        
        fprintf([' Aggregate Markup at ', num2str(round(Mu,2)), '\n'])
        
        fprintf('\n');
        
        fprintf(' A. Assigned Parameter \n');  
        
        fprintf('\n');
        
        disp(table(symbol, description, round(values, 2), 'VariableNames', {'Symbol','Variable Description','Value'}))
        
        fprintf('\n');

        fprintf(' B. Calibrated Parameter \n');
        
        fprintf('\n');
        
        symbol      = categorical({'M';'-';'-';'b';sprintf('%c ', 958);sprintf('%c ', 963);[sprintf('%c', 949),'/',sprintf('%c ', 963)];sprintf('%c ', 981)});
        description = categorical({'aggregate markup';'top 5% sales share';'materials share';'regression coefficient';'Pareto tail';'demand elasticity';'super-elasticity';'weight on valued-added'});
        values      = [Mu; moment_model(4); B; moment_model(6); p.xi; p.sigma; moment_model(6); p.phi];
        
        disp(table(symbol, description, round(values, 2), 'VariableNames', {'Symbol','Variable Description','Value'}))
        
        fprintf('\n');

    case 'T3'

        fprintf('<strong> Table 3: Markup Disperion and Productivity Losses, Benchmark Model </strong> \n'); 

        fprintf('\n');
        
        fprintf([' Aggregate Markup at ', num2str(round(Mu,2)), '\n'])

        fprintf('\n');
        
        description = categorical({'25th percentile';'50th percentile';'75th percentile';'90th percentile';'99th percentile'});
        values      = [wprctile(mu, 25, weight); wprctile(mu, 50, weight); wprctile(mu, 75, weight); wprctile(mu, 90, weight); wprctile(mu, 99, weight)];
        
        fprintf(' Cost-weighted distribution of markups: \n');  
        
        fprintf('\n');
        
        disp(table(description, round(values, 2), 'VariableNames', {'Variable Description','Value'}))
        
        fprintf('\n');
        
        description = categorical({'gross output ';'value added';'value added, M = 1'});
        values      = [-(Z/Zp - 1)*100; -(GDP/GDPp  - 1)*100; -(GDP1/GDPp - 1)*100];
        
        fprintf(' Aggregate productivity loss in percent: \n');  
        
        fprintf('\n');
        
        disp(table(description, round(values, 2), 'VariableNames', {'Variable Description','Value'}))
        
        fprintf('\n');

    case 'T4'

        switch experiment

            case {'uniform', 'entry'}

                start_subsidy

            case {'efficient', 'sizedep'}

                start_planner

        end

end


